
> # FILE:     UnmaskingReplication.R
> # AUTHOR:   Austin Hart & J. Scott Matthews
> # DATE:     July 7, 2021
> # SUBJECT:  REPLICATION FILES, JOURNAL OF POLITICS
> #           Main Analysis for Exps 1 - 5
> 
> # Notes for replication:
> #   Required packages: tidyverse, knitr, broom
> #   User-written functions & graphical theme defined below
> 
> 
> # SETUP ---------------------------------------------------
> 
> # load packages
>   library(tidyverse)
-- Attaching packages ------------------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.3     v purrr   0.3.4
v tibble  3.0.6     v dplyr   1.0.4
v tidyr   1.1.2     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.1
-- Conflicts ---------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

>   library(knitr);library(broom)

> # directory (set directory below if you choose)
>   # setwd("your.path.here")
> 
> # Data for main text, experiments 1-5    
>   load(url("https://www.dropbox.com/s/c9safb3z4bib6wb/HartMatthewsJOP.RData?dl=1"))

> # CUSTOM FUNCTIONS ---------------------------------------
>   
> # FctClean() - function to relabel + reorder factors 
>   FctClean <- function(facvar, ...){
+     fct_recode(fct_relevel(facvar, ...), ...)
+   }  

> # Custom theme for visuals
>   mytheme = theme_bw() + theme(
+     plot.title = element_text(hjust = 0, size = 10, face = 'bold'),
+     legend.title = element_text(size = 10),
+     legend.text = element_text(size = 10),
+     axis.text.x = element_text(size = 10),
+     axis.title.x = element_text(margin = unit(c(3, 0, 0, 0), "mm"), size = 10),
+     axis.text.y = element_text(size = 10), axis.title.y = element_text(margin = unit(c(0, 3, 0, 0), "mm"), size = 8),
+     panel.grid.minor = element_blank(),
+     panel.grid.major.y = element_line(color = 'gray80',linetype = 'solid',size = 0.35),
+     panel.grid.major.x = element_blank(),
+     axis.ticks = element_blank(),
+     panel.spacing = unit(2, "lines"),
+     plot.margin=unit(c(.2,.2,.2,.2),"cm"),
+     panel.border = element_blank(),
+     strip.background = element_blank(),
+     strip.text = element_text(size = 10, face = 'bold')
+   )    

> # RETENTION FUNCTIONS -------------------------------------
> 
> # Filter for Table 4
>   df =
+     dfmain %>%
+     select(VoteB,AvgPre,AvgPost,Study,Design)

> # Table 4: estimated performance vote by study
>   df %>%
+     split( f = list(.$Study,.$Design), drop = TRUE ) %>%
+     map_df(
+       ~ broom::tidy( lm(VoteB ~ AvgPre + AvgPost, data = .x) ),
+       .id = 'Study'
+     ) %>%
+     select( -std.error, -statistic ) %>%
+     filter( term != '(Intercept)' ) %>%
+     pivot_wider(
+       values_from = c( estimate, p.value ),
+       names_from = c( term )
+     ) %>%
+     mutate_at( # scale performance in 100s of units
+       vars(starts_with( 'estimate' )),
+       ~ . * 100
+     ) %>%
+     mutate_at( # 1-sided tests
+       vars(starts_with( 'p.' )),
+       ~ . * 0.5
+     ) %>%
+     arrange(Study) %>%
+     .[, c(1, 2, 4, 3, 5)] # reorder columns
# A tibble: 7 x 5
  Study         estimate_AvgPre p.value_AvgPre estimate_AvgPost p.value_AvgPost
  <chr>                   <dbl>          <dbl>            <dbl>           <dbl>
1 1.Unmasking           -0.0929       3.74e-19            0.135        1.18e-21
2 2.Recognition         -0.0554       2.15e- 6            0.123        3.83e-19
3 3.Unmasking           -0.114        9.09e- 7            0.158        6.98e- 8
4 4.Recognition         -0.0650       2.77e- 6            0.153        5.16e-22
5 4.Unmasking           -0.129        1.89e-20            0.151        4.30e-16
6 5.Recognition         -0.0466       3.63e- 5            0.156        1.04e-37
7 5.Unmasking           -0.138        4.06e-26            0.192        5.09e-29

> # summary: baseline retention rate
>   df %>% # by study
+     group_by( Study ) %>%
+     summarise( Reappointed = 100 * mean(VoteB) )
# A tibble: 5 x 2
  Study Reappointed
* <dbl>       <dbl>
1     1        75.9
2     2        71.5
3     3        73.8
4     4        68.1
5     5        64.1

>   mean(100 * df$VoteB) # overall
[1] 69.50557

> # remove objects  
>   rm(df)  

> # STUDY 1 ---------------------------------------
>   
> # Isolate data
>   df1 = 
+     dfmain %>%
+     filter(Study == 1)

> # Table 5: Unmasking incumbent performance
>   # Estimates shown
>     df1 %>%
+       lm(VoteB ~ 0 + NoiseS1 + I(AvgAB/100):NoiseS1 + I(AvgA/100):NoiseS1, data = .) %>%
+       summary()

Call:
lm(formula = VoteB ~ 0 + NoiseS1 + I(AvgAB/100):NoiseS1 + I(AvgA/100):NoiseS1, 
    data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0268 -0.1143  0.1517  0.2773  0.5980 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
NoiseS1High(250)               0.22264    0.20060   1.110    0.267    
NoiseS1Low (100)               0.27804    0.21036   1.322    0.187    
NoiseS1High(250):I(AvgAB/100)  0.11653    0.01828   6.374 3.03e-10 ***
NoiseS1Low (100):I(AvgAB/100)  0.16414    0.02127   7.716 3.38e-14 ***
NoiseS1High(250):I(AvgA/100)  -0.07117    0.01328  -5.358 1.08e-07 ***
NoiseS1Low (100):I(AvgA/100)  -0.12500    0.01600  -7.815 1.64e-14 ***
---
Signif. codes:  0 �***� 0.001 �**� 0.01 �*� 0.05 �.� 0.1 � � 1

Residual standard error: 0.401 on 843 degrees of freedom
Multiple R-squared:  0.7895,	Adjusted R-squared:  0.788 
F-statistic: 526.9 on 6 and 843 DF,  p-value: < 2.2e-16


>   # Tests: differences/interaction effect 
>     df1 %>%
+       mutate(NoiseS1 = relevel(as.factor(NoiseS1), ref = 'Low (100)')) %>%
+       lm(VoteB ~ I(AvgAB/100)*NoiseS1 + I(AvgA/100)*NoiseS1, data = .) %>%
+       summary() # divide p-values by 2 for 1-tailed tests    

Call:
lm(formula = VoteB ~ I(AvgAB/100) * NoiseS1 + I(AvgA/100) * NoiseS1, 
    data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0268 -0.1143  0.1517  0.2773  0.5980 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    0.27804    0.21036   1.322  0.18662    
I(AvgAB/100)                   0.16414    0.02127   7.716 3.38e-14 ***
NoiseS1High(250)              -0.05540    0.29068  -0.191  0.84888    
I(AvgA/100)                   -0.12500    0.01600  -7.815 1.64e-14 ***
I(AvgAB/100):NoiseS1High(250) -0.04761    0.02805  -1.698  0.08996 .  
NoiseS1High(250):I(AvgA/100)   0.05383    0.02079   2.589  0.00979 ** 
---
Signif. codes:  0 �***� 0.001 �**� 0.01 �*� 0.05 �.� 0.1 � � 1

Residual standard error: 0.401 on 843 degrees of freedom
Multiple R-squared:  0.1281,	Adjusted R-squared:  0.123 
F-statistic: 24.78 on 5 and 843 DF,  p-value: < 2.2e-16


>   rm(df1)  

> # STUDY 4 ---------------------------------------
>   
> # Isolate Study 4 data
>   df4 =
+     dfmain %>%
+     filter(Study == 4) %>%
+     select(Design,Report1.S4,Report2.S4)

> # Figure 1: a preference for benchmarking
>   # Name and mutate
>     df4 =
+       df4 %>%
+       mutate_all(as.character) %>%
+       rename(
+         Rep1 = Report1.S4,
+         Rep2 = Report2.S4
+       )

>   # Tabulate report selections 
>     iso =
+       df4 %>%
+       pivot_longer(
+         -Design,
+         names_to = 'name',
+         values_to = 'Choice'
+       ) %>%
+       group_by(Design,Choice) %>%
+       summarise(n = n()) %>%
+       mutate(
+         prop = 100 * n/sum(n),
+         grp = 'Individual selections'
+       )
`summarise()` has grouped output by 'Design'. You can override using the `.groups` argument.

>   # Identify & tabulate choice combinations    
>     com =
+       df4 %>%
+       mutate(
+         Choice = case_when(
+           Rep1 == 'B' & Rep2 != 'AvB' ~ 'Per',
+           c(Rep1 %in% 'A' & Rep2 %in% 'B') | c(Rep1 %in% 'AvB' | Rep2 %in% 'AvB') ~ 'Ben',
+           TRUE ~ 'Irr'
+         )
+       ) %>%
+       group_by(Design,Choice) %>%
+       summarize(n = n()) %>%
+       mutate(
+         prop = 100 * n/sum(n),
+         grp = 'Combination'
+       ) %>%
+       ungroup()
`summarise()` has grouped output by 'Design'. You can override using the `.groups` argument.

>   # Append and factor
>     com =
+       com %>%
+       bind_rows(iso) 

>     com$Choice = as_factor(com$Choice)

>     com =
+       com %>%
+       mutate(
+         Choice = FctClean(
+           Choice,
+           'Incumbent avg' = 'B',
+           'Comparator avg' = 'A',
+           'Contrast' = 'AvB',
+           'Bonus to date' = 'Y',
+           '(none)' = 'None',
+           'Inc-centered' = 'Per',
+           'Benchmarked' = 'Ben',
+           'Irrelevant' = 'Irr'
+         )
+       )

>   # Visualize
>     f1 = 
+       com %>%
+       ggplot( aes(x = Choice, y = prop )) +
+       facet_grid( ~ fct_rev(grp), scales = 'free', space = 'free' ) +
+       geom_point(aes(shape = fct_rev(Design)), 
+                  color = 'black', size = 2.75, position = position_dodge(width = 0.25)) +
+       labs(
+         x = NULL,
+         y = 'Percent who chose the report',
+         color = 'Task design',
+         shape = 'Task design'
+       ) +
+       scale_y_continuous(expand = c(0,0),
+                          breaks = seq(0, 100, 50),
+                          labels = c('0','50','100%')) +
+       coord_cartesian(ylim = c(0,103)) +
+       scale_shape_manual(values = c(24,19)) +
+       mytheme +
+       theme(
+         axis.text.x = element_text(angle = 90,vjust = 0.4,hjust=1)
+       )

>   # Export  
>     ggsave(filename = "fg1.png", plot = f1, 
+            width = 6, height = 3, dpi = 1200)  

> # model estimates by report chosen
>   # Isoalte Data and Vars from main
>     df4 = 
+       dfmain %>%
+       select(Design,Report1.S4,Report2.S4,VoteB,AvgPre,AvgPost) %>%
+       mutate_at(c('Report1.S4','Report2.S4'),as.character)

>     df4 =
+       df4 %>%
+       mutate(
+         Choice = case_when(
+           Report1.S4 == 'B' & Report2.S4 != 'AvB' ~ 'Per',
+           c(Report1.S4 %in% 'A' & Report2.S4 %in% 'B') | c(Report1.S4 %in% 'AvB' | Report2.S4 %in% 'AvB') ~ 'Ben',
+           TRUE ~ 'Irr'
+         ),
+         choiceComp = if_else(Choice == 'Ben','Benchmarker','Else')
+       )

>     df4 %>%
+       split(list(.$Design,.$choiceComp)) %>%
+       map_df(
+         ~ tidy(lm(VoteB ~ AvgPre + AvgPost, data = .x)),
+         .id = 'Arm.Choice'
+       ) %>%
+       select(-std.error,-statistic) %>%
+       filter(term != '(Intercept)') %>%
+       pivot_wider(
+         values_from = c(estimate,p.value),
+         names_from = c(term)
+       ) %>%
+       mutate_at(
+         vars(starts_with('estimate')),
+         ~ 100 * .
+       ) %>%
+       arrange(Arm.Choice) %>%
+       .[, c(1, 2, 4, 3, 5)] # reorder columns  
# A tibble: 4 x 5
  Arm.Choice              estimate_AvgPre p.value_AvgPre estimate_AvgPost p.value_AvgPost
  <chr>                             <dbl>          <dbl>            <dbl>           <dbl>
1 Recognition.Benchmarker         -0.0909       1.08e- 6            0.157        8.77e-15
2 Recognition.Else                -0.0487       6.44e-10            0.143        1.73e-61
3 Unmasking.Benchmarker           -0.123        7.22e-17            0.143        4.47e-13
4 Unmasking.Else                  -0.114        3.98e-50            0.161        1.34e-56

> # clean enviro   
>   rm(f1,com,iso,df4)

> # STUDY 5 ---------------------------------------
>   
> # Isolate data
>   df5 = 
+     dfmain %>%
+     filter(Study == 5) %>%
+     select(Study,Design,VoteB,AvgPost,AvgPre,ReportS5)

> # Figure 2: can targeted cues mitigate benchmarking?  
>   # Estimate models
>     me = # extraction game
+       df5 %>%
+       filter(Design == 'Unmasking') %>%
+       lm(VoteB ~ 0 + ReportS5 + I(AvgPost/100):ReportS5 + I(AvgPre/100):ReportS5, data = .) 

>     mr = # recognition game
+       df5 %>%
+       filter(Design == 'Recognition') %>%
+       lm(VoteB ~ 0 + ReportS5 + I(AvgPost/100):ReportS5 + I(AvgPre/100):ReportS5, data = .) 

>   # Store estiamtes/ses
>     emods =
+       tibble(
+         game = c(rep('Unmasking',6),rep('Recognition',6)),
+         report = rep(c('Incumbent Avg','Pre-Benchmarked','Bonus (control)'),4),
+         vars = c(rep('Incumbent Tenure',3),rep('Pre-Incumbent',3),
+                  rep('Incumbent Tenure',3),rep('Pre-Incumbent',3)),
+         coefs = c(me$coefficients[4:9],mr$coefficients[4:9]),
+         ses = c(sqrt(diag(vcov(me)))[4:9],sqrt(diag(vcov(mr)))[4:9])
+       ) %>%
+       mutate(
+         game = fct_rev(game),
+         report = factor(report, 
+                         levels = c('Incumbent Avg','Pre-Benchmarked','Bonus (control)'))
+       )

>   # Plot
>     f2 =
+       emods %>%
+       ggplot(aes(x = fct_rev(report), y = coefs, 
+                  fill = fct_rev(vars), shape = fct_rev(vars))) +
+       facet_grid(game ~., scales = 'free_y') +
+       geom_linerange(aes(ymin = coefs - 1.96*ses, ymax = coefs + 1.96*ses),
+                      size = 1) +
+       geom_point(size = 3, color = 'black') +
+       coord_flip() +
+       geom_hline(yintercept = 0) +
+       mytheme +
+       scale_fill_manual(values = c('white','black')) +
+       scale_shape_manual(values = c(21,22,24)) +
+       labs(x = NULL,
+            y = 'Effect of performance\non Incumbent retention',
+            color = NULL,
+            shape = NULL,
+            fill = NULL) +
+       theme(panel.grid.major.y = element_blank(),
+             panel.grid.major.x = element_line(color = 'gray80'),
+             legend.position = 'right')

>     ggsave(filename = 'fg2.png', plot = f2, height = 3, width = 5.5, dpi=1200)

>   # Estimated performance vote by design & by treatment arm     
>     df5 %>%
+       mutate(Report = relevel(ReportS5, ref = 'Bonus to date')) %>%
+       lm(VoteB ~ I(AvgPost/100)*ReportS5 + I(AvgPre/100)*ReportS5, data = .) %>%
+       summary()

Call:
lm(formula = VoteB ~ I(AvgPost/100) * ReportS5 + I(AvgPre/100) * 
    ReportS5, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0293 -0.4287  0.1690  0.3453  0.8888 

Coefficients:
                                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          -0.234819   0.209242  -1.122   0.2620    
I(AvgPost/100)                        0.147753   0.015674   9.426  < 2e-16 ***
ReportS5Benchmark                     0.202826   0.292737   0.693   0.4885    
ReportS5Bonus to date                 0.194553   0.285269   0.682   0.4954    
I(AvgPre/100)                        -0.071604   0.013851  -5.169  2.7e-07 ***
I(AvgPost/100):ReportS5Benchmark      0.024183   0.022675   1.066   0.2864    
I(AvgPost/100):ReportS5Bonus to date  0.001079   0.022381   0.048   0.9616    
ReportS5Benchmark:I(AvgPre/100)      -0.042877   0.019794  -2.166   0.0305 *  
ReportS5Bonus to date:I(AvgPre/100)  -0.024162   0.019762  -1.223   0.2217    
---
Signif. codes:  0 �***� 0.001 �**� 0.01 �*� 0.05 �.� 0.1 � � 1

Residual standard error: 0.432 on 1367 degrees of freedom
Multiple R-squared:  0.1945,	Adjusted R-squared:  0.1898 
F-statistic: 41.25 on 8 and 1367 DF,  p-value: < 2.2e-16


>     df5 %>%
+       mutate(Report = relevel(ReportS5, ref = 'Bonus to date')) %>%
+       lm(VoteB ~ I(AvgPost/100)*ReportS5 + I(AvgPre/100)*ReportS5, data = .) %>%
+       summary()

Call:
lm(formula = VoteB ~ I(AvgPost/100) * ReportS5 + I(AvgPre/100) * 
    ReportS5, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0293 -0.4287  0.1690  0.3453  0.8888 

Coefficients:
                                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          -0.234819   0.209242  -1.122   0.2620    
I(AvgPost/100)                        0.147753   0.015674   9.426  < 2e-16 ***
ReportS5Benchmark                     0.202826   0.292737   0.693   0.4885    
ReportS5Bonus to date                 0.194553   0.285269   0.682   0.4954    
I(AvgPre/100)                        -0.071604   0.013851  -5.169  2.7e-07 ***
I(AvgPost/100):ReportS5Benchmark      0.024183   0.022675   1.066   0.2864    
I(AvgPost/100):ReportS5Bonus to date  0.001079   0.022381   0.048   0.9616    
ReportS5Benchmark:I(AvgPre/100)      -0.042877   0.019794  -2.166   0.0305 *  
ReportS5Bonus to date:I(AvgPre/100)  -0.024162   0.019762  -1.223   0.2217    
---
Signif. codes:  0 �***� 0.001 �**� 0.01 �*� 0.05 �.� 0.1 � � 1

Residual standard error: 0.432 on 1367 degrees of freedom
Multiple R-squared:  0.1945,	Adjusted R-squared:  0.1898 
F-statistic: 41.25 on 8 and 1367 DF,  p-value: < 2.2e-16


> # clean objects
>   rm(f2,emods,me,mr,df5)
